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ABSTRACT 

The interaction rates of dark-matter halos and subhalos, including collisions and 
mergers, are computed using high-resolution cosmological iV-body simulations of the 
ACDM model. Although the number fraction of subhalos of mass > 2 x 10 11 /i -1 Mq 
is only ~ 10%, we find that the interaction rate of such subhalos is relatively high 
because they reside in high density environments. At low redshift, the subhalo collisions 
dominate the total collision rate, and even at z = 3 they are involved in more than 
30% of all collisions. About 40% of the "major" mergers (those of mass ratio > 0.3) 
are between subhalos. Therefore subhalo interactions must be incorporated in models 
of structure formation. 

We find that the collision rate between halos in physical density units, is oc (l+z) s , 
with 5 = 3 — 4, in agreement with earlier simulations and most observational data. 

We test previous analytic estimates of the interaction rates of subhalos based 
on statistical models, which could be very inaccurate because of the small number 
of subhalos and the variation of conditions within small host halos. We find that, 
while such statistical estimates may severely overestimate the rate within hosts < 
10 13 /i _1 M Q , typical of high redshifts, th ey are valid for large r hosts regardless of the 
number of subhalos in them. We find the Making & Hut (1997) estimate of the subhalo 
merger rate to be valid for hosts <; 6 x lO il /i~ 1 Af at all redshifts. 

The collision rate between subhalos and the central object of their host halo is 
approximated relatively well using the timescale for dynamical friction in circular 
orbits. This approximation fails in ~ 40% of the cases, partly because of deviations 
from the assumption of circular orbits (especially at low redshift) and partly because 
of the invalidity of the assumption that the subhalo mass is negligible (especially at 
high redshift). 



1 INTRODUCTION 



Interactions among dark matter halos and the galaxies that 
reside i n them are inevitable in any cosmological scenario 
of hierE rchical structure formation and have important ob- 



Galaxy interactions are observed directly (e.g., Davi. 
ct al. 1997j |dc Mcllo et al. 1997j ; [Conselice fc Gallagher 199" 



Ncuschaefer et al. 1997| ) and indirectly by statistical means 
and t hrough their outcomes (|Abraham, Ellis, fc Glazebrook 



1997 



servati qnal implications^ JHjy^Mrrteract^^ tidal strip- 

ping ar d harassment ( Mocirej^ake^fc^ They 



prompt ; the exchange of angular momentum (e.g., 
fc Efsta thiou 1987 ), trigger star bursts (e.g., | Mihps 



Barnes 



Hern 



quist 1994a| ; |Mihos fc Hernquist 1994b ), mix different stel- 



l ar populations (Larson fc Tinsley 197£), produce tidal tails 
£ roo mre fc Toomre 1972|; pubinski, Mihos, fc Hernquist 



inter a ction rate rises as fun c tion of redsh i ft (|Le Ffevre et al 
200(1; patton fc et al. 200C|; patton 1998j parlberg, Couch 



man, fc Thomas 199ty |Zcpf fc Koo 198S| ) (but cf. parTbcTg 



Conselice, Bershady, fc Jangrcn 2000| ; |Wu 1999|). The 



et 



al. 200C| ) and the density of the environment. Since the in- 



199£ 



S^rmg^T^^^hite_l_99_9[ ) , and change the galaxy mor- 



phology (e.g., Demler 1992). The tidal stripping and ha- 



rassment that occur during these interactions enrich the in 



terga lactic medium with processed elements (Kolatt et al 



200C), disrupt existing galaxies, build new ones, and heat 
up the galaxy interior matter. 



teraction rate is roughly proportional to the number density 
squared times the relative velocity, groups and clusters of 
galaxies show a higher specific interaction rate (per object) 
than the specific rate in the field. 

The details of the physical processes that take place 
during and as a result of these interactions are not well un- 
derstood. Gas dynamics, star formation and supernova feed- 
back hinder a direct understanding of how such interactions 
influence the process of galaxy formation in general. In or- 
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der to overcome this difficulty, semi-anal y tic investigations 
dKauffmaim, White, & Guidcrdoni 1993j |Cole et al. 1994 



Somerville & Primack 1999; Somerville, Primack, & Faber 



200C ) u 3e simplified models for the merger rate of substruc- 
ture (including the role of dynamical friction) in order to 
explore the possible implications of these events. These sim- 
plified models should be tested and improved by comparison 



with re ilistic cosmological simulations. 



But even before addressing the additional physics in- 
volved in galaxy interactions, the halo gravitational interac- 



tion rale alone is too complicated to calculate analytically. 



The finite halo size and density profile, the many body in- 
teractions, the varying background potential all have their 
impact on the interaction rate. 

The first stage of the theoretical exploration of gravi- 
tational halo interact ions introduces the impulsive approxi- 
mation (Bpitzer 1955), to estimate the effect of a close pas- 
sage of two gravitationally bound systems. This approxima- 
tion models the heating up of each individual system by 
the transfer of gravitational potential energy to internal ki- 
netic energy. It also models the disruption by stripping and 
the contraction of the remain ing debris. The impulsive ap- 



proxim ation has been tested (Richstone 1975 



Dekel, Lecar 



& Shalam 1980) in iV-body simulations where the effects 
of close encounters on mass exchange and profile evolution 
were addressed. The main objective of these investigations 
was to test mass loss due to tidal stripping. The conclu- 
sion was that although the impulsive approximation works 
well even beyond "fast" encounters, it is invalid for encoun- 
ters with small impac t parameters. Later, usi ng two-halo 
iV-body experiments, Augilar fc White (1985 ) found that 
the impulsive approximation works only if the minimal dis- 
tance between the passing halos is bigger than ~ 10 effective 
(de Vaucouleurs) radii. This clearly prevents the use of this 
approximation for close encounters and collisions. 

The evolutio n of many-gal axy systems (rich clusters) 
was addressed by Merritt (1983) with the inclusion of close 



encounters and tidal stripping, but not merging. Merritt in- 
vestigated the question of mass segregation between subha- 
los and the host cluster, and found that the implied dynam- 
ical friction time is longer than would be expected without 
the presence of close encounters and tidal stripping. 

Merging of dark m atter halos was investi gated on the 
basis of simulations by |Makino & Hut (1997|) (MH), who 



parameterized equal-mass mergers as a function of impact 
parameter and relative velocity for different initial density 
profiles. MH also explored the parameter space of very close 
passages. 

Recently, many authors have performed detaile d simu- 



et al. 1996; Tormen 1997; fTormen, Bouchet, & White 19971 


Tormen 1998; Ghigna et al. 1998: sail 


den Bosch et al. 


199£ 
199£ 




ensui, Funato, & Makino 199E ; |Okamoto & Habe 
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bore et al. 1999(). These investigations focused on is- 



sues like tidal stripping, harassment, and orbital parameters. 
They led to various insights about the role interactions play 
within clusters, understanding how morphology changes due 
to fast encounters, how anisotropic orbits affect the final 
density profile, quantifying the dynamical equilibrium state 
in clusters, and quantifying the mass fraction of subhalos 
with respect to the diffuse matter. However, systems the 
size of rich and poor galaxy groups have been neglected un- 



til now. The present work is the first detailed investigation 
of subhalo encounters within a cosmological volume simula- 
tion. 

Previously, the only relevant cosmological modeling of 
halo interactions has been via ana lytic methods such as th e 
Press-Schechter (PS) formalism (|Press fc Schcchtcr 1974|) 



and its exten ded versions (EPS) (Lacey & Cole 1993; Lacey 



Cole 1994). These theoretical predictions have been com- 



pared to iV-body results, wit h surprisingly good agreement 



(up to a fa ctor of ~ 2) (e.g., Gross et al. 1998 



Somerville 



et al. 200C). Unfortunately EPS theory only accounts for 
mergers between distinct, virialized halos, and does not con- 
sider any other form of interaction. In particular EPS does 
not take into account non-merger interactions and does not 
provide the joint multiplicity-mas s function of the progeni- 
tors (cf. Somerville fc Kolatt 19991) . In addition, EPS theory 
neglects the spatial distribution of the existing halos, en- 
vironmental dependence, and, most importantly, the inter- 
actions between subhalos within a given encompassing viri- 
alized host halo. All of these cases are extremely relevant 
to galaxy formation and evolution. For example, the results 
of MH were used by Somerville et al. (2000) to include in 
semi-analytic models the process of subhalo mergers, which 
they considered to be related to the L yman break galaxie s 
observed at high redshift (z <; 3). In Kolatt et al. (199!:) 
we emphasize the role these collisions play, without appeal- 
ing to MH but rather as identified in the high-resolution 
simulations studied in more detail here. 

One may think that the shortcomings of the PS/EPS 
formalism may be overcome by using EPS to model isolated 
systems and then applying kinetic theory to virialized sys- 
tems and substructure (cf. Somerville et al. 2000) . However, 
this approach is questionable - first because of the cross- 
talk between isolated halos and subhalos, and then because 
of the limitations of the kinetic theory. 

In the kinetic theory within the framework of statisti- 
cal mechanics one predicts the collision rate among many 
bodies in a given system. When we calculate collision or 
merger rates in a gravitating many body system, the simple 
approach of kinetic theory is no longer valid. Several impor- 
tant differences prevent us from applying the conventional 
methods of statistical mechanics: 



(i) The gravitating systems under consideration are not 
necessarily in equilibrium. The velocity distribution function 
may be anisotropic, may deviate from the usually assumed 
Maxwellian distribution for relaxed systems, and may evolve 
in time. Moreover the number of "bodies" (i.e., halos) is not 
constant; halos may merge or dissolve and other new halos 
may enter the system or just pass through it in the course 
of time. 

(ii) The "bodies" are not equal, but exhibit a large range 
of masses. The mass function (and thus the cross section 
distribution) is complicated and evolves in time and space. 

(iii) The "bodies" are not point-like or rigid bodies, nor 
do they retain their identity; "bodies" may shrink, evapo- 
rate, grow, deform, etc. 

(iv) Gravitating system s exhibit "gravitational focusing" 
(Binney & Tremaine 1987), i.e., an increase of the cross sec- 
tion due to the gravitational attraction between the inter- 
acting bodies. 

(v) Most of the gravitational systems under consideration 
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are embedded within a background potential (e.g., individ- 
ual galaxy halos in groups or clusters) . This potential cannot 
in general be neglected. Moreover, the interplay betwee n the 
cluster or group potential and its subhalos (see item (iii) ) 
evolves in time. 

(vi) The effect of the background potential sometimes 
adds anisotropy to the orbits - for example spiraling-in of 
s atellite halos within a ho st halo due to dynamical friction 
(Binney & Tremaine 1987), or radial orbits along filaments 



and to the center of deep potential wells. 

For the reasons outlined above, at present the only way 
to measure the halo interaction rate and other character- 
istics is to resort to high-resolution iV-body simulations. 
There are some minimum requirements for any simulation 
used to tackle such a problem. One needs a big enough dy- 
namical range in spatial and mass resolution for the com- 
pilation of a statistical sample of the host/subhalo systems 
within a cosmological volume. The mass ratio between sub- 
halos and hosts can range from ~ 1 to three or four or- 
ders of magnitude. The evolution of number density of the 
host halos and their subhalos is strongly dependent on the 
cosmological background and its ability to provide a fresh 
supply of progenitors throughout the host evolution. In ad- 
dition, spatial resolution is of paramount importance since 
substructure tends to erase within the host unless the reso- 
lution allows its clear identification. Another necessary con- 
dition is the storage and analysis of many time-steps of the 
simulation to determine the temporal behavior and interac- 
tion of the substructure. 



The si mulations we use in the current analysis (Klypin 



et al. 1999) meet all of these requirements, and in doing so 



represent an important advance in simulation technology. 
We ha ve designed, te s ted, and used a new hierarchical halo 



finder (Bullock 1999; Bullock et al. 200C hereafter BKS+ 



specifically for the goals of this paper. 

In this paper we focus on the DM component alone and 
analyze the interactions between halos and within host ha- 
los - those that have at least one other halo within their 
virial radius. We make no discrimination by host or sub- 
halo masses, though we sometimes segregate the results into 
clusters and groups. A major goal of this paper is to assess 
various estimates for the time scales in virialized systems 
that consist of distinct constituents. 

In section ^| we briefly discuss the iV-body simulations 
and describe our strategy to find and model the dark mat- 
ter halos. We shall then define the hierarchy we build for 
halos within the simulation and the various criteria for in- 
teractions. Section ^ prepares the ground for the interaction 
rate calculation by discussing the multiplicity function of 
host halos and its sensitivity to mass resolution. The mul- 
tiplicity function plays an important role in various analyt- 
ical estimates for the collision rate. In we shall present 
the collision rate of halos, and the progenitor mass spec- 
trum. A comparison to various approximations for collision 
timescales is carried out in and we devote ^ to a simi- 
lar evaluation of the dynamical friction timescale. We finally 
summarize our findings in §(jj. 



2 SIMULATED HALOS 

We used the Adaptive Refinemen t Tree (ART) code 



(Kravtsov, Klypin, & Khokhlov 1997) to simulate the evo- 
lution of collisionless DM in the "standard" ACDM model 
(fi m = 1 - fi A = 0.3; H = WOh = 70 km s" 1 Mpc" 1 ; 
as = 1.0). The simulation followed the trajectories of 256 3 
particles within a cosmological periodic box of size L = 
60 /i -1 Mpc from redshift z = 40 to the present. A basic 
512 3 uniform grid was used, and up to six refinement levels 
were introduced in the regions of highest density, implying 
a dynamic range of ~ 32, 000. The formal resolution of the 
simulation is thus / rcs ~ 2h _1 kpc, and the mass per DM 
particle is m p ~ 1 x 10 9 hT 1 Mq. We analyze 15 saved out- 
puts at times between z — 5 and z — 0. 

A complementary, L = 30/i _1 Mpc simulation of eight 
times higher mass resolution and two times higher force reso- 
lution is used to check issues of completeness and robustness 
to resolution. We analyze ~ 10 time-steps of this simulation 
in the range 1.7 < z < 7. 



2.1 Model halos 

The identification of halos is a key feature of the analysis; We 
try to make it objective and self-con siste nt when following 
halo interactions and mergers (see §2.4). Traditional halo 



finders utilize either friends-of-friends algorithms or over- 
densities in spheres or ellipsoids to identify virialized halos. 
These algorithms fail to identify substructure. We there- 
fore have designed a new hierarchical halo finder, based 



on the bou nd density maxima (BDM) algorithm (Klypin 
et al. 19991). The detail s of the halo finder are described 
elsewhere (|Bullock 1999, BKS+), and we summarize below 



only its main relevant features. 

We impose a minimum number of 50 particles per mod- 
eled halo, unify overlapping maxima, and iteratively find the 
centers of mass of spheres about the maxima. We compute a 
spherical density profile about each center and identify the 
halo virial radius R v i r inside which the mean overdensity has 
dropped to a value A v j r , based on the spherical infall model. 
For the family of flat cosmologies ( fl m + Qa = 1), the value 
of A v i r can be approximated by (Bryan & Norman 199£ ) 
A vir ~ (18vr 2 + 82a; - 39z 2 )/(l + x), where x = (l(z) - 1. 
In the ACDM model used in the current paper, A v i r varies 
from about 200 at z > 1 to A v i r — 340 at z = 0. If an 
upturn occurs in the density profile inside -Rvir, we define 
there a truncation radius Rt- 

An important step of our procedure is the fit of the den- 
sity profile out to the radius min(_R v i r , R t ) with a universal 



functional form . We adopt the NFW profile (Navarro, Frenk 
& White 19951), 



PNFw(r) = 



Pa 



(r/R s )(l + r/R a y 



(1) 



with the two free parameters R s and p 3 — a characteris- 
tic scale radius and a characteristic density. This pair of 
parameters could be equivalently replaced by other pairs, 
such as Rs and -R v i r - The minimum halo mass correspond- 
ing to 50 particles is ~ 5 x W 10 h~ 1 M Q (for the 60 h~ 1 Mpc 
simulation). The modeling of the halos with a given func- 
tional form allows us to assign to them characteristics such 
as a virial mass and radius, to estimate sensible errors for 
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Figure 1. Various relationships among halos and their nomen- 
clature. The thin lines represent i? v i r and the shaded regions - 



these quantities, and to remove unbound particles. Using 
the NFW fits, we iteratively remove unbound particles from 
each modeled halo and unify every two halos that overlap 
in their i? s and are gravitationally bound. Finally, we look 
for virialized regions within R a of big halos to enable sub- 
halos near the centers of big host halos, e.g., mimicking cD 
galaxies in clusters. 

Our halo finding is complete only for halos contain- 
ing more than ~ 200 particles, i.e., of mass larger than 
~ 2 x 10 11 h~ 1 M e (for the 60 7i -1 Mpc simulation) or ~ 
1.5 x 10 10 A _1 M Q (for the 30 ft _1 Mpc simulation). The re- 
sulting halo catalog is therefore gradually incom plete below 
these masses, as described in Bullock (1999) and |5igad et al 



(2000|Qhereafter SKB+). This should be kept in mind be- 
low, when we sometimes refer also to masses somewhat be- 
low the full completeness limit. 



2.2 Hierarchy classification 

Two halos are classified as a subhalo and a host when the 
center of the smaller halo lies within i? v i r of the larger one. 
The one exception is the case we term "partners" , in which 
the two centers lie within R v i r of each other but the small 
halo is not fully contained within the large one and the mass 
ratio is smaller than 4/3. 

In the current application we limit ourselves to the first 
level of subhalos (but the classification scheme can straight- 
forwardly be extended to deal with many levels of subhalos 
within subhalos and with other relations between halos). 
All halos that are not subhalos are termed "distinct" halos; 
they can be either "isolated" halos or hosts. Halos that do 
not host any subhalos are termed "simple" halos; these can 
be either isolated halos or subhalos. Figure |l| illustrates the 
different possible relationships and their classification. 



2.3 Interaction definitions 

The definition of halo interaction in an iV-body simulation 
is somewhat arbitrary. We can easily identify violent inter- 
actions, but more subtle phenomena are less evident. The 
borderline between a merger and a close encounter, for in- 
stance, should depend upon the final bound/unbound stage 
of a two body system, but finding a final bound state at some 
later timestep does not necessarily imply that a merger had 
occurred. Since in reality, as in the simulations, halos can 
disintegrate, or lose some of their mass during an interac- 
tion, the decision of what we would call a merger, namely 
what fraction of the progenitor halo mass should merge with 
another progenitor halo to form a united halo, is a matter 
of definition. 

In our procedure we followed these guidelines: (a) All 
the definitions should be reproduceable within any N- 
body simulation that is being analyzed, (b) Avoid arbitrari- 
ness as much as possible in the interaction definitions, (c) 
Avoid sensitivity of the results to particular choices of in- 
teraction definition. The decision to model all halos by one 
unique functional form was motivated by these guidelines. 



2.4 Finding a collision 

Interactions, by their nature, are a dynamical process. Any 
definition of interaction must therefore involve the compar- 
ison of two or more time-steps in the simulation. We hence 
identify and model halos in stored time-steps of the sim- 
ulations, and compare the simulated halos in subsequent 
time-steps in order to define interactions. 

Our procedure for finding a collision that occurred dur- 
ing an interval Az = 22 — zi (22 > 21) is as follows. We list 
all particles that are bound to the halo hi at 22. Particles 
are only listed under the highest level halo they are bound 
to - i.e., if under a subhalo, then not also its host halo. We 
identify these same particles at the later time step z%, if they 
are listed under any halo at this stage. 

If any of the identified halos at z\ contains particles 
from hi we fit an NFW profile based on these particles only. 
The fitting parameters are listed under the halo hk(zi). 

If at the end of this process a halo hk(zi) has more than 
one set of fit parameters from previous redshift halos, we 
check whether any pair of sets obey the "collision criterion" 
(see below) . Note that halos that at 22 contained more than 
twice the minimal number of particles, can merge with two 
separate halos at a later timestep 21 < 22, although this is 
unlikely to happen. 

The collision criterion is consistently defined to match 
the halo-finder; if two density profiles, as calculated for two 
sets of particles, overlap in their scale radius _R S , a collision 
has occurred. If each of the two modeled sets is gravitation- 
ally bound to the other modeled set then the collision is a 
merger. If this is not the case, the collision is an "unbound" 
collision. We use this term and the term "free-free collision" 
interchangeably. Since each fit contains an error in the de- 
rived Rs value, the occurrence of a collision is up to the 
combined error in _R S of the two fits. The binding relation is 
defined by comparing the relative center of mass velocity to 
the modeled NFW escape velocity of hi at the location of 
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3 STATISTICAL VIEW ON SUBSTRUCTURE 



Most galaxies in the local universe (J> 50%) are in groups ac- 
cording to various observational definitions of groups (e.g., 
Zabludoff & Mulchaey 1998 and references therein) . If these 
definitions coincide with physically bound virialized sys- 
tems, and if one identifies the simulated dark matter halos 
with galaxies (up to a certain mass) then we should find a 
similar fraction of subhalos in the simulations, provided the 
lower mass cutoff of the simulations accommodates the same 
galaxy mass range for which the group analyses have been 
carried out (see SKB+). 

In our standard simulation we identify a total of ~ 7300 
halos at z = and ~ 4500 at z = 4. Between 5% to 15% 
of all halos in the simulations we use are subhalos, almost 
independent of cosmological epoch and the mass resolution 
of the simulation. According to our classification, about 10% 
of all subhalos are subhalos of level 3 and higher. 

This leads to ~ 15 — 20% of so-called "galaxies" be- 
ing group members (including the host if it is not too mas- 
sive) which is less than the observed value for the simulation 
mass range. The mismatch between the identified fraction 
of grouped galaxies versus the same for grouped halos is re- 
solved (SKB+) by allowing a halo tree buildup in which one 
uses three times the halos' virial radius to identify subha- 
los, namely replace i? v j r by 3/? v ir in the definitions (i)-(iv) in 
52.2. However in the current context where we are interested 



l°MM h0 , t /(h-iM s )] > 13 



in real virialized systems, we use the notion of substructure 
defined in 



2:2 



One of the important ingredients for the collision prob- 
ability function is the number density of subhalos within 
a given halo. This number depends first and foremost on 
the mass of the host halo, but the cosmological epoch, the 
power spectrum, and the host halo shape may also affect it. 
At first we focus on the host mass. Figure [] shows the multi- 
plicity function for all host halos with Mhost > 10 13 /i -1 -Mq 
at z = 0. We can identify three large clusters in the sim- 
ulation (TVsub = 17, 25, 29 at z — 0) but we cannot detect 
whether the multiplicity function is continuous throughout 
the subhalo number range due to insufficient statistics. The 
multiplicity function is a function of the halo/subhalo mass 
resolution. Using the same simulation w ith > 2 times better 
halo mass resolution, |Colin et al. (1999 ) found 4 to 10 times 
more subhalos for the same host halos. Note that Colin et al. 
are defining halos differently and do not need to construct 
a density profile for the halos in order to identify collisions. 
When a mass cutoff similar to the one used here is applied to 
the higher mass resolution halo catalog, a very similar multi- 
plicity function is obtained. For an effective mass cutoff ten 
times smaller than the 100% completeness value used here 
(2 x 10 11 /i -1 Mq) the multiplicity function for hosts more 
massive than 10 13 /i _1 Mq peaks at 4-5. As we shall shortly 
see, in all of the following comparisons with model predic- 
tions we shall use the actual halos we identify and hence 
we avoid the need to use only halos that have masses larger 
than our 100% completeness limit. The completeness func- 
tion for our halo finder is presented elsewhere (SKB+). We 
are ~ 70% complete for the smallest mass halo used in this 
analysis. 

For all host halos in the current analysis, almost regard- 
less of their mass, the multiplicity function peaks sharply at 
the first bin, namely at a situation of a binary system host 



E . 

Z 10 - 



MM^Ah-'Mj] > 13 
z = 2.9 



Figure 2. The multiplicity function of all host halos with mass 
> 10 13 /i _1 M Q at z = (top) and z = 2.9 (bottom). Note the 
different axis scale. Recall that complete identification is obtained 
only for halo masses of M > 2 X 10 11 ?t- 1 MQ. 



and a single subhalo. For comparison in Fig. we plot 
the same multiplicity function at z = 2.9 where there are 
no prominent clusters present, and the multiplicity declines 
sharper than the z = function for iV su b = 2 — 4. The mul- 
tiplicity function does not seem to be a very strong function 
of the mass resolution of the simulation for high mass hosts. 
A somewhat more detailed view of the multiplicity-mass de- 
pendence is shown in Fig. ^| where a scatter plot reveals the 
general expected trend of high multiplicity for more massive 
halos. It is clearly seen though that on occasion even rela- 
tively low mass halos (< 3 x 10 13 /i _1 Mq) exhibit high mul- 
tiplicities. Figures ^a and |^b show the multiplicity function 
of somewhat less massive hosts at z = 2.9 for the standard 
resolution and for eight times better mass resolution, respec- 
tively. We lowered the host mass threshold to 1O 12 /i _1 M 
at both resolutions in order to improve on the statistics. 
The limiting factor is the small number of massive halos in 
the small volume of the high resolution simulation rather 
than the mass resolution itself. This problem becomes more 
severe at higher redshifts. 
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MM^tAh-'Mj] > 12 
z = 2.9 
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Figure 3. Scatter plot in the multiplicity 
z = host halos. 



host mass plane for 



Figure |B| shows the average number of subhalos as a 
function of redshift for 2 mass bins, M less than or greater 
than 1O 13 /i _1 M0. The dashed lines show the average num- 
ber for all the halos, regardless of whether they contain sub- 
halos or not, whereas the solid lines focus only on halos that 
contain at least one subhalo. Most halos in the larger mass 
bin contain subhalos. The average number of subhalos for 
host halos rises substantially from ~ 1.0 to ~ 2.0 between 
z = 5 and z = 2 and levels to a constant of ~ 2.3 subhalos 
per host for z ~ — 1. In the small mass bins the common 
situation is one subhalo per host, i.e., a binary system. 



l°g[M h0 , t /(h-'M s )] > 12 
z=2.9 

High mass resolution 



4 COLLISION RATE 

Recall that halos that overlap in their NFW R 3 parameters, 
and that mergers are bound collisions (§^J). In addition 
to determining the collision rate of all halos, in this paper 
we concentrate on the collision rate among subhalos, where 
the multiplicity function dictates much of the characteris- 
tic behavior. We may already guess from the lack of strong 
redshift evolution of the structure of this function (Fig. bk 
in comparison to Fig. ^b and Fig. ^ that the subhalo colli- 
sion rate does not decrease monotonically with time at late 
times, as expected for the collision rate among distinct halos 
in an e xpanding universe. Figure p|a shows the collision rate 



Figure 4. The multiplicity function of all host halos with mass 
> lO 12 ^" 1 Mq at 2 = 2.9 for the lower resolution simulation 
(top) and for the higher resolution (bottom), uncorrected for the 
different volumes. 



resolution effects). Only a fraction of these collisions occur 
among subhalos or between hosts and their subhalos. The 
shaded area on Fig. |^a depicts the collision rate only for col- 
lisions that involve at least one subhalo. A frequently quoted 
statistic is the merger rate (what we call the collision rate) 



in physical units as a fu nction of redshift ( Carlberg, Couch 



as a Au ction of look-back time or equivalently the redshift. man, fc Thomas 1990 ; Zepf & Koo 198E; Patton & et al 



ume is described by the solid line. The errors are derived 
from the errors in the modeled Rs of the colliding halos 
t hro ugh the definition of a collision as overlapping 7? s (cf. 
§2.4). Note that the collision and merger rates with finer 



The overall number of collisions per Gyr per comoving vol- 200C ; Lc Fcvrc et al. 200C ). Figure [j] shows this representa- 



time resolution were calculated, and in every case found to 
be consistent with the values that appear in Fig. ^ (and 
the following figures) within the quoted errors. In the over- 
lapping redshift range of the standard and high resolution 
simulation, given the same mass cuts we find similar rates 
within the errors. That reduces the uncertainty with regard 
to the effect of the mass resolution on the calculated rate 
for a given mass range (see below for a discussion of other 



tion of the collision rate where the oc (1 + z) increase at 
low redshift is prominent. Similar increases are reported by 
the abovementioned observational analyses in the z = — 1 
range, although not by Carlberg et al. (2000) - suggesting 
that selection effects and the different observational tests for 
collisions used by different authors may be important. 

In Fig. ^| we zoom-in onto the shaded region of Fig. |^a. 
Collisions that involve subhalos can further be divided into 
subclasses. The shaded region in Fig. ^ shows the merger 
(as opposed to all collisions) rate among subhalos, while the 
solid thick line describes the host - subhalo collision rate, 
namely a subhalo that got "absorbed" by its host or went 
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Figure 5. The average number of subhalos per host (triangles, 
solid lines) and per halo (circles, dashed lines) in two mass bins: 
log[M/(h~ 1 JWq)] < 13.0 and > 13.0 (lower line and upper line 
respectively). Top: standard resolution, bottom: high mass reso- 
lution (eight times standard). 



through its center without keeping its identity. If the subhalo 
keeps its identity it is classified as a "cD" and therefore does 
not show up on the collision budget. A subhalo, however, can 
merge with a "cD" that resides at the center of a host. 

As expected the overall collision rate rises steeply from 
z ~ 4 to z ~ 2, reflecting the picture of bottom up hierar- 
chical clustering and merger history. Below z = 1 the rate 
starts to drop gently because in a low-f2 m ACDM cosmol- 
ogy the accelerating expansion halts the collisions/mergers 
of isolated virialized halos. However in this redshift range 
the collisions among subhalos start to become a substantial 
fraction of all collisions. By that time galaxy groups and 
clusters start to form, and a significant increase in the num- 
ber density of subhalos is obtained in these big host halos. 
The collision rate increases correspondingly. The increased 
velocity dispersion helps the process as well, though free- free 
collisions are still less probable than mergers. The fraction 
of unbound collisions out of all collisions among subhalos re- 
mains almost constant. Most of the erratic temporal features 



Figure 6. The collision rate of dark matter halos with (top) M > 
5 X lO lo ft _1 M (complete only above 2 X 10 11 /i~ 1 M Q ) (bottom) 
M > 6 X lO 9 fc- 1 M (complete only above 2.5 X 10 10 /i _1 Af a ) 
as function of look-back time. All collisions are described by the 
solid line where the error bars are model errors (see text). The 
fraction of collisions that involve at least one subhalo is marked 
by the shaded area. 



in the collision rate are due to unbound collisions between 
a subhalo and its host. 

The collision rate is sensitive to the mass resolution of 
the simulation and its completeness limit. The higher mass 
resolution simulation not only shows a much higher collision 
rate (Fig. ^b) as a result of the mass function steepness, but 
it also shifts the peak of the collision rate to higher redshift 
where small halos are already virialized and participate in 
the collisions either as isolated halos or as subhalos of hosts. 
Figure ^b shows the overall collision rate and subhalos col- 
lision rate at a resolution eight times higher than the one in 
Fig. ^ only down to z = 1.7. At z = 2 the overall collision 
rate with the high mass resolution is ~ 7 times higher and 
the rate among subhalos ~ 10 times higher. Recall that the 
product halo mass function reaches its completeness 1 value 
only at ~ 200 particles, i.e., 2 x 10 11 ft -1 M for the stan- 
dard resolution. This problem is less severe than it could be 
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Figure 7. Same as figure |£j, here in physical units which empha- 
size the increase of the collision rate and as a function of log(l+z). 
The two reference lines are 7^(0)(l + z) 3 and 7£(0)(1 + z) 4 , where 
TZ(0) is the rate at z = 0. 
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Figure 8. The collision rate of dark matter subhalos as a func- 
tion of look- back time (thin solid line). The total number corre- 
sponds to the shaded area in Fig. where the shaded area here 
is the merger rate among subhalos. The thick solid line is the 
collision rate between subhalos and hosts. 



because the slope of the low mass end of the subhalo mass 
function (SKB+) is only ~ —0.7 to —0.5. The mass spec- 
trum of the progenitors in the two resolutions is of course 
very different. For any specific comparisons with observa- 
tions (e.g., Kolatt et al. 1999) one has to first figure out the 
relevant mass range for progenitors. 



Collisions 0.03>z>0.00 
log[M h /(h-'M )]> 10.73 
Complete @ 11.3 

~ Distinct halos 
* Subhalos 



Bigger progenitor log(M/h 'M s ) 



Figure 9. Progenitor mass spectrum for collisions between z = 
0.03 and 2 = 0. Empty circles mark masses of progenitors that 
are both distinct halos and solid triangles represent progenitors 
of collisions that involve at least one subhalo. Here and in the 
next figures, these masses are based on the particles that actually 
participate in the collision. 



4.1 Progenitor mass spectrum 

Two important issues have to be examined with respect to 
the progenitor mass distribution. It is important to learn the 
mass range of progenitors in order to relate them eventually 
to physical entities. In particular, whether a collision is a 
"major" or "minor" one depends on whether the progeni- 
tor mass ratio is greater or less than some particular value, 
typically 0.3. 

Figure [] shows the distribution of masses for pairs of 
progenitors of collisions. The quoted masses of the progeni- 
tors are calculated directly from the number of particles that 
took part in the collision, which may be only fractions of the 
progenitor halo, as opposed to the modeled masses (of either 
the progenitor halos or the product halo). These masses can 
therefore be smaller than the mass resolution limit which 
refers to the minimal halo mass. In most cases these masses 
are a very good representation of the noisier modeled pro- 
genitor mass, especially for large mass progenitor halos. 

As shown in Fig. |^a, the fraction of collisions among 
subhalos at low redshift is high. In Fig. ^ we show the mass 
spectrum for all the collisions between 0.03 > z > 0.00 
divided into collisions between distinct halos (circles) and 
among subhalos (triangles). At this low redshift almost all 
bigger progenitors have already exceeded ~ 10 11 /i _1 Mq. 
While most of the progenitors of collisions among distinct 
halos have mass ratios of ~ 1, the collisions among subha- 
los involve a much broader spectrum of mass ratios. The 
complexity of the distribution clearly demonstrates why the 
numerical approach is essential. The mass function of sub- 
halos cannot be derived by EPS (cf. SKB+), nor can the 
progenitor mass function. Most triangles at the right end of 
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Figure 10. Progenitor mass spectrum for collisions between z ~ 
4 and z ~ 3. Symbols are as in Fig. ^. (a) The low resolution 
simulation (b) High resolution simulation. 



the figure are a subhalo - host collision. The pile-up of pro- 
genitors at the low mass end is a reflection of the steepness 
of the mass function (cf. Somerville & Kolatt 1999; SKB+). 
Figure [lO^ , is similar to Fig. |^ but at a higher redshift 



(3.9 < % < z.y). At this higher redshift the collisions among 



subhalos bt»ai a much greater resemblance lu Lhe piugeni- 

tor mass spectrum of the distinct halos. The higher reso- 
lution simulation (Fig. |l"c|b) at similar redshifts shows the 
general trend of smaller mass progenitors which were not 
identified in the standard resolution simulation, and similar 
spectrum coverage for distinct halo collisions and collisions 
among subhalos. There is also a redshift evolution (SKB+) 



of the 



tlubhalb rnass function that folds Into the prbgeriitor 
mass fuiicliun aL differ enL ledshifL 



4.2 Major vs. minor mergers 

The mass ratio between merger progenitors is assumed 



in semi-analytic models (Kaiifhiiaim, White, fc Guider- 
doni 19931; bole et al. 1994 ISomerville fc Primack 



Somerville, Primack, fc Fabcr 200C ) to determine the mor- 



19991 
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Figure 11. Probability distribution of mass ratio for progenitors 
of mergers among subhalos (top) as a function of the redshift 
of the mergers. The redshift range is indicated on the right. 
The same distribution for all mergers (bottom). The probability 
distributions for all collisions are very similar. 



phology of the central galaxy within the merger product. 
Base d on hydrodynamic simulations (Barnes & Hernquist 



1992), it is assumed that the central galaxy resulting from a 
"major merger" (dark halo mass ratio greater than 0.3) will 
have a dominant spheroidal stellar component, while the re- 
sult of a "minor merger" of a small galaxy (mass ratio less 
than 0.3) with a larger disk galaxy is again a disk galaxy. If 
indeed this identification is valid, we expect a higher frac- 
tion of "major mergers" in a high density environ ment, for 
consistency with the density-morphology relation (Dressier 



198C). It is therefore interesting to see whether the mass 
ratio distribution alters when we turn to examine only col- 
lisions among subhalos. 

Figure |Tl] a shows the probability distribution of the pro- 
genitor mass ratio for mergers among subhalos, at different 
epochs. At high redshift (lower part) the distribution is al- 
most uniform, and toward lower redshifts it tends to con- 
centrate at low mass ratios, namely accretion of small pro- 
genitors onto bigger subhalos. The extended set of all colli- 
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sions among subhalos shows a similar behavior. If mergers 
among distinct halos are also included (Fig. [TTJa) there ap- 
pears at lower redshifts another probability peak for similar 
mass progenitors. 

If we take mass ratios 0.3 or greater as leading to major 
mergers, we find that ~ 50 — 60% of all mergers are major. 
Mergers among subhalos start with being major in ~ 60% of 
all cases at z ~ 4 and only ~ 40% by z = 0. Similar results 
are obtained if unbound collisions are included. 



5 TESTING INTERACTION RATE 
APPROXIMATIONS 

In section [I] we have listed the various reasons for the inade- 
quacy of the traditional analytic calculation for the collision 
rate as derived in statistical mechanics. In this section we 
shall attempt to compare results from the statistical me- 
chanics approach and from two-halo interaction simulations 
with the full iV-body results obtained here. 



5.1 Analytical calculations 

In the kinetic theory, the average number of collisions a par- 
ticle experiences per unit time is the collision rate, 1Z co ii- 
For a system in equilibrium TZ co u = naV, with n the 
number density of the particles, a their collisional cross- 
section, and V the average relative velocity between two 
particles in the ensemble. For a Maxwellian velocity distri- 
bution V = Wrms/V^j where « rms is the one-dimensional 
velocity dispersion. If the particles interact gravitationally, 
and are all identical in mass (to) and physical size (r C oll), 
the effective cross-section increases (gravitational focusing, 
e.g., Binney & Tremaine 1987) and the collision rate changes 
correspondingly by the addition of 16 G rnn r co u /V with 
r 'coii — o'/tt. For a heterogeneous system with particles that 
retain their spherical symmetry, but exhibit various r co n and 
to values, we therefore define the average collision rate by 

— 16 Gmr co u 1 



{R-coii) =n[aV-\ = J 



(2) 



These rates are calculated for rigid bodies with radii r co H> 
and yield a result for a collision rate as opposed to the 
merger rate. 

Makino & Hut (1997, MH) derived an analytic approx- 
imation for the merger rate, based on their binary collision 
simulations. MH tested for mergers of equal-mass isolated 
halos of various initial density profiles. Instead of perform- 
ing a big N-body simulation, they covered a large portion of 
phase space with their simulations and derived their expres- 
sion from the results. As MH pointed out, there are numer- 
ous limitations with this approach. The merging criterion in 
MH is the requirement of negative energy value when the 
two systems are treated as a two body system. There is no 
simple way to determine particle affiliation with halos after 
the collision occurred (though the MH procedure seems very 
reasonable) . 

MH concluded that there is only a weak dependence on 
the initial halo profile and that a good approximation for 
the merger rate in a multi-halo system would be 

\ ( -Rl/2sub X 



n 



2 x 10~ 3 iV 2 



-Rho 

1 Mpc J \ 0.1 Mpc 



(^rms sub \ / 
100 kms-V V 



^rms host 

300 kms- 1 



Gyr" 



(3) 



This approximation should be relevant for a closed system 
(host halo) of radius i?h os t which contains N identical subha- 
los with half mass radius i?i/2sub each. The one-dimensional 
velocity dispersion for the host and the subhalo are v lms host 
and frmssub, respectively. Note that MH neglected in this 
approximation the likelihood that the central object of their 
"cluster" may have a substantially larger cross section. Here 
we include it among the subhalo population of the host. 
Obviously the MH approximation is only applicable to our 
work if generalized to a system with appropriate progenitor 
mass spectrum and varying cross sections. 



5.2 Comparison with simulations 

When we calculate the various statistical predictions, using 
the data of the simulations and the formulae of the analytic 
or semi-analytic approximations, we are bound to introduce 
an error. While the approximations explicitly assume no evo- 
lution for the subhalo population and neglect the host poten- 
tial, the number density and density profiles of the subhalos 
actually change from one time step to another, and the back- 
ground potential of the host changes too. In comparison to 
realistic iV-body results a compromise should be taken. We 
would like the time lag between two compared timesteps to 
be as short as possible in order to validate the no-evolution 
assumption and the neglect of second order effects. On the 
other hand, we would like a long enough time-step to en- 
able accumulation of proper statistics. The number density 
of subhalos is not constant because of the entry of new sub- 
halos in between the analyzed time-steps, and subhalo de- 
struction due to tidal stripping and galaxy harassment. 

We here check how well these approximations compare 
with the results obtained from the simulations. Namely, we 
compare the modified kinetic theory approximation with our 
measured collision rate, and the MH approximation with the 
measured merger rate. At each time step we record for each 
host halo the number N of subhalos it contains, and calcu- 
late the number density n = iV/VW using the modeled virial 
volume of the host. For consistency with our collision defi- 
nition (cf. §2.4), we take the cross section for each subhalo 
to be a function of its modeled R s . The instantaneous rate 
for this host halo is evaluated via 



E 



n (R s ) V + ■ 



V 



w 



(4) 



where the sum is over all subhalos. W % = RI/o-r, s are the 
weights for the modeled 7? s 's of the subhalo population 
((W 1 ) 2 for the simple naV estimate). The predicted number 
of collisions for this host during the proper time At between 
subsequent registered redshifts is obtained by 



N cM = N s TZ co u At , 



(5) 



with N 3 the number of subhalos of the host halo. 

For statistical reasons it is more robust to compare the 
global rate rather than the rate for each individual host halo. 
Moreover, halos do not keep their identity from one time- 
step to another. 

Figure [l2] shows the average number of mergers and col- 
lisions per host halo per Gyr as a function of look-back time. 
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Figure 12. The average number of (a) mergers and (b) collisions per host halo per Gyr as a function of look-back time. The errors for 
collisions are bigger due to larger analysis errors on the modeled R B of the free-free collisions. 



The temporal evolution of the collision rate per host shows 
a similar behavior to that of the overall rate per comoving 
volume element (Fig. |^) with an increasing rate until z ~ 2 
(z ~ 1 for mergers) a decline until z — 0.2 and another rise 
thereafter. Note that the comoving merger rate per host halo 
remains close to constant throughout the redshift evolution 
and only varies in the range 0.5 — 2Gyr -1 . The same mild 
variation applies to the collision rate per host halo, which 
remains in the range 1 — 2Gyr _1 . 

Figure [llj shows the same data as appear in Fig. [l^ 
in comparison to the approximation predictions. The empty 
symbols of Fig. jL3| represent the approximations for the col- 
lision (merger) rate per host according to Eqs. ^-|^. Note 
that the input data are taken from the true parameters of 
the simulation halos (Eq. ||). 

Similarly, taking into account the fact that we are deal- 
ing with a realistic mass spectrum, and not equal-mass ha- 
los, and conforming to the half mass radius for the NFW 
modeled profiles, we use the same procedure (Eqs. ^, |^) to 
calculate the expected rate a la MH for each host halo in 
the sample. We summed all rates (multiplied by the number 
of subhalos for each host) and obtained the MH prediction 
for the merger rate per host halo. 

ft is evident from inspection of Fig. |l^ that the theoret- 
ical approximations do not agree with the measured quan- 
tities at high redshift, and that the disagreement tends to 
decrease at lower redshifts. Part of the reason for the dis- 
crepancy can be read directly from Fig. ^| — even massive 
halos tend to have only a handful of subhalos sufficiently 
massive to survive our cuts, and do not reach the required 
large number density for statistical considerations. Other 
reasons for the discrepancy will be discussed shortly. Note 
that gravitational focusing only slightly changes (by 0—50%) 
the straightforward predictions of kinetic theory. 



A way to cope with the limited number of subhalos per 
host is to divide the host population by the host mass or the 
number of subhalos they host. Although the approximations 
break down, they may still work for a subset of halos where 
the approximation assumptions are valid. In the next section 
we explore this possibility. 

5.3 Interactions as a function of iVs U b and host 
mass 

Wc first test the conjecture that the governing factor for the 
mismatch between the predictions of the approximations is 
the small number of subhalos. At each given time-step the 
available multiplicity function does not allow, in general, a 
large enough number of host halos with high multiplicity. 
In order to check the conjecture, though, one needs proper 
statistics, i.e., many host halos of high multiplicity. A way 
around this difficulty is found by stacking redshift steps to- 
gether according to the multiplicity of the host halo under 
consideration. That allows us to accumulate statistics from 
all time-steps and bin the data by iV su b of the halos instead 
of their collision redshift. The error we introduce by this 
statistics accumulation is mainly due to different conditional 
mass function of subhalos at different redshifts (SKB+). 

Figure |l^ shows the average number of mergers (col- 
lisions) per host halo per Gyr as a function of the num- 
ber of their subhalos. We compare only host halos with 
M > 2 x 10 11 since at high redshift both the host halo and 
the subhalo of lower masses have very low concentration (cf. 
Bullock et al. 2000) , and thus may obey the "cD" case crite- 
rion and yield unphysically high rates. Evidently the rate is 
still overestimated with all three approximations by about 
one to two orders of magnitude. A fixed logarithmic shift 
of the kinetic theory model by ~ —1.6 — 1.7 for N su h ^ 25 
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Figure 13. The average number of (a) mergers and (b) collisions per host halo as a function of look-back time and redshift, in comparison 
with analytic models. Solid symbols describe results from the current analysis (errors smaller than dot size). Empty circles on the left 
describe the MH predictions of merger rates within host halos, empty symbols on the right are obtained from the collision rate prediction 
of kinetic theory (triangles) and the additional gravitational focusing (squares). 



provides reasonably good agreement between the measured 
collision rate and the calculated one, for the high resolution 
simulation as well. We were able to identify the cause for the 
increasing discrepancy at the high multiplicity bin. It comes 
from low mass hosts with high multiplicity at Z2 (cf. Fig. [}]) 
and some higher mass hosts with high multiplicity at Z2 . For 
the high mass hosts that show the most discrepant result (a 
small fraction), the multiplicity changes dramatically from 
one time-step to the next, generally toward lower multiplic- 
ity at the later time. The approximations (Q) and (5) take 
N s , the number of subhalos, and n, their number density, to 
be constant. But in these discrepant cases N 3 decreases due 
to tidal destruction, causing the subhalos to drop below our 
resolution limit, and thus n decreases as well. 

In conjunction with the additional increase in R v i T , the 
approximations show much higher interaction rates than are 
detected in the simulations. The breakdown of the approxi- 
mation assumptions suggests that there might yet be a dif- 
ferent factor which governs the collision rate more strongly 
and causes the deviation from the predictions. 

The other governing parameter for the approximation 
success may be the host mass. Massive host halos on aver- 
age have higher multiplicity and thus better fulfill the sta- 
tistical assumptions. In addition to that, massive hosts are 
more stable and their gravitational potential prevent subha- 
los from escaping the virial radius. The competing effect is 
the stronger tidal force they exhibit; however, this is some- 
what relaxed due to their typical lower concentration param- 
eter (cf. BKS+) and hence moderate density gradients. We 
use the same stacking approach by which we examined the 
approximation dependence on the multiplicity to circumvent 
the difficulty of not having enough statistics at each time 
step. Here we consider mass bins, accumulating all similar 
host masses from all time steps. Since the halos were cho- 
sen under the proviso they had already virialized, similar 



host masses correspond to similar velocity dispersion within 
them, with no regard to the cosmological epoch, although 
this is not strictly valid with respect to the number of sub- 
halos at different epochs and their scale radii (BKS+). 

Figure [j^ shows the merger rate as a function of the 
host mass. The simulation statistics only allow 2 — 3 signif- 
icant mass bins. Clearly for the high mass bin the MH ap- 
proximation works well, while it overestimates the rate by a 
factor ~ 5 for the low mass bin, for reasons discussed in the 
previous subsection. The collision rate prediction from the 
kinetic theory works well for the large mass bin as well and 
exhibits a discrepancy of about an order of magnitude for 
the low mass bin. These results explain the large deviation 
from the models seen, e.g., in Fig. |l3| The mass function is 
always dominated by the low mass end of the hosts. Natu- 
rally, this end tends to consist of less massive hosts at higher 
redshift. Most hosts therefore do not fall in the regime where 
good agreement with the models is obtained, and the aver- 
age rate per host is dictated by these. The conclusion is that 
the difference between the kinetic theory predictions with 
or without taking gravitational focusing into account is not 
big, and they both fit the results from the simulations rather 
well for masses larger than ~ 10 13 h -1 Mq. In the high reso- 
lution simulation a similar trend is detected and successful 
prediction of all three models can be extended due to better 
statistics to M host Z 1O 12 ' 5 /i _1 M . 

In an attempt to bin the mass more finely a useful ap- 
proximation for all masses of host halos is found to be 

simulation ; (6) 

with log(/) = -0.27 log(M/10 18 ' 5 /i _1 M Q ). The approxima- 
tion (El) works fine for both mass resolutions when applied 
to the kinetic theory (naV, with or without gravitational 
focusing) . There is no similar correction formula for the MH 
predictions. 
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Figure 14. The average number of (a) mergers and (b) collisions 
per host halo as a function of iV su b (the number of subhalos of the 
host) for host halos > 2 X 10 11 h^ 1 Mq, in comparison to analytic 
models. Symbols are as in Fig. |l3j 



Figure 15. The average number of (a) mergers and (b) collisions 
per host halo as a function of the host mass in comparison to 
analytic models. Symbols are as in Fig. |l^. The lower "gravita- 
tional focusing" point on the higher mass bin reflects the different 
weighting scheme. 



In spite of the poorer statistics, we now go back to 
the evolution of the collision rate in time (Fig. |l3| ) and ap- 
ply the same calculation under the constraint of Mhoat > 
6.3 x 10 11 hT 1 Mq. This allows statistically significant results 
at almost all redshifts. Figure shows the results. A very 
nice agreement is obtained for all redshift bins between the 
MH approximation and the simulation. The kinetic theory 
still overestimates the collision rate by about an order of 
magnitude beyond redshift z — 1, but gets closer to the 
simulation results at lower redshifts if gravitational focusing 
is ignored. 

The MH prediction is expected to give somewhat higher 
results than the simulation since MH calculations consider 
the half mass radius while in the simulation the requirement 
is of a cross section R s in conjunction with a final bound 
state. Figs. |l5| a and p"6| a show that on average for large host 
halos the two merger definitions are similar, but in details 
and for small halos they differ. It is encouraging to verify 
that the merger definition is quite robust and the outcome 
rate does not depend on its fine details. The kinetic theory 



predictions tend to overestimate the rate since they neglect 
stripping (and thus shrinking cross section) and harassment 
to the degree of disappearance of halos. 



6 DYNAMICAL FRICTION 

The breakdown of the kinetic theory statistics in the limit of 
a few-body system, and in particular for low masses as seen 
in the simulations, leads us to seek an alternative assess- 
ment for the collision rate in these systems. Fig. |l7] shows 
the distance of the identified collision, i? C oiiision from the 
host center, for all collisions that involve at least one sub- 
halo. Note that beyond the first radial bin, the rise in the 
number of collisions as a function of the radius in which 
they occur is slower than expected for a homogeneous dis- 
tribution throughout the host's virial volume, as assumed 
in the kinetic theory approximations (Eq. ^, |^). Instead, the 
distribution follows more closely a mass profile of M tx r. 
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Figure 16. Same as figure [1^, except that here only host halos of mass > 6.8 X 10 11 /! 1 Mq are included in the comparison to obtain 
a much better match with the Makino & Hut approximation. Note that for the kinetic theory the discrepancy becomes smaller too. 



In the limit of a binary system, every collision is by defini- 
tion a host - subhalo collision. These are the host - subhalo 
collisions that contribute to the big spike in Fig. [l7| near 

-Rcollision/^vir = 0. 

Two competing processes are responsible for this type 
of collisions: a ra dial orbit may lead the subh alo directly to 
the host's center (Moore, Lake, & Katz 1998), while a more 



tangential orbit may end up spiraling into the center of the 
host due to dynamical friction. The characteristic time for 
spiraling in, as estimated by dynamical-friction theory, is 
clearly relevant for the latter but it may or may not be a 
useful approximation for the former. We therefore set to test 
the applicability of the dynamical-friction time scale in the 
realistic systems of our simulation. 

Traditionally the rate of collisions is estimated using a 



simplifi ed form of dynamical friction ( Kauffmann, White 



Guiderdoni 1993 



Somerville fc Primack 1999| ) which uti- 



lizes a very simple approach. The halo is assumed to be an 
isothermal sphere, and the subhalo is assumed to be a point 
mass on circular orbits, ignoring mass -loss as it spirals in . 
The more elaborate approximations by Benson ct al. (200C) 
and by Klypin et al. (1999) allow a more realistic density 
profile for the halos but they still assume circular orbits for 
subhalos. 

We here focus on collisions between subhalos and their 
hosts. The progenitor subhalo of these collisions is identified 
in the regime between R V [ T and R s of the host (apart from 
the "cD" case). We first check how well the traditional dy- 
namical friction timescale predicts the behavior in the simu- 
lation. The expression for the dynamical friction time in the 
fra mework of NFW density p rofile for the host halo is given 
by (iBinney & Tremaine 19871) 



NFW 
Tdf 



tdtjr) ( d In Mhoat (r) 



2r 



dlnr 



+ 1 dr 



(7) 



ix q r 




Figure 17. Distribution of the collision locations in units of the 
host iJ v j r . The diagonal line describes the relative volume in the 
radial bins, normalized by the counts at the second inner-most 
bin. 



where ( Klypin et al. 1999 ) 



47rG 2 (lnA)M sub p(r) [erf(X) - 2X e - x2 /^\ 

vl{r) = GM hoat (r)/r; X = v c /(V2a r ); 
i?virM host (r) 



A = 



rM sub 



(8) 
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and oy is the one-dimension velocity dispersion at radius r, 



SIS 



2 2 2x(l+xf 

/(x) = ln(l + a;)- 

Wmax = « c (2.15iis) . 



x 3 (l + x) 2 

T 



dx, 



1 + x ' X R s 



(9) 

The halo profile, and therefore the circular velocity 
v c (r), are assumed to remain constant between subsequent 
time-steps. The radii n and rj are the initial and final radii 
to consider, respectively. The initial separation is identi- 
fied at the earlier redshift as the offset between the center 
of mass of the host and that of the subhalo. The final radius 
Tj is the sum of the R s values for the host and the subhalo, 
to maintain consistency with the collision definition. For the 
"cD" case 77 is taken to be zero. The choice for the subhalo 
mass is somewhat ambiguous: the scheme we use modeles 
the subhalo mass according to its NFW density profil e. A n 
alternative is to consider only the mass within R t (c.f. §2.1). 
We tried both possibilities and found only negligible differ- 
ences between the two results. Note, however, that the sub- 
halo mass does change in between the two time-steps and 
thus deviations from the no-mass-loss approximation (Eq. 
0) are expected. Interestingly, since the subhalo mass usu- 
ally decreases, the approximation should underestimate the 
dynamical friction time. We shall soon see that the opposite 
happens, in many cases, indicating that the simplified model 
may be crude for other reasons too. 

An alternative approach would be to calculate raf in a 
similar way it is done in certain semi-analytic models (e.g., 
Somerville & Primack 1999). The host halos are assumed 
to follow a singular isothermal sphere (SIS) profile whose 
constant circ ular velocity is given by the virial mass ad ra- 
dius: v c = \/GM v ir/ R v i T . The dynamical frictio n time scale 
is then simplified to (Binney & Tremaine 1987) 



sis 



1.17 (r? 



rj) v c 



In A GAf sub ' (10) 
where the Coulomb logarithm is approximated by A = 
M hoBt /M Buh . 

Among all collisions that have taken place between Z2 
and Z\ we select only collisions of the type subhalo - host, 
namely systems that at Z2 already belonged to the same viri- 
alized system where approximations for dynamical friction 
time are applicable. We identify the progenitor that at Z2 
used to be the subhalo and calculate its expected tos accord- 
ing to Eqs. (Q) and ([to]). We do not take its actual orbital 
parameters into account, because we are interested in find- 
ing out whether r<jf is a good approximation when applied 
to a statistical sample of collisions as in the semi-analytic 
models. We then compare the expected dynamical friction 
time of this individual collision with the elapsed time At 
between Z2 and z\ in the simulation. The time resolution of 
the analyzed simulation output times does not allow a de- 
tailed comparison of the calculated dynamical friction time 
and the collisions as they happened in the simulations. The 
proper time gap between successive stored times serves in- 
stead as an upper limit for the real time. 

If the dynamical friction time is a reliable estimate, all 
the individual r^t values should be smaller than At, because 
we selected for all collisions that have actually happened on 
a time shorter than or equal to At. The choice of going 




-4 -2 



Log(r d[ /At) 



NFW 




Log(T dr /At) 

Figure 18. Distribution of dynamical friction timescale T<jf in 
units of the elapsed time At between subsequent simulation out- 
puts for various redshifts (right). Here r,jf is calculated assum- 
ing that the host halo radial profile is either singular isothermal 
sphere (top panel) or NFW (bottom panel). 



"backwards" in time, namely from identified collisions to 
their progenitors, is made because some of the subhalos at 
the higher redshift Z2 got disrupted by the time of the lower 
redshift Z\ and can no longer serve as dynamical friction 
probes. 

Figure [l8| shows the distribution of dynamical friction 
time in units of At, for all progenitor subhalos in the simu- 



© 0000 RAS, MNRAS 000, 000-000 



16 Kolatt et al. 



lation at a given redshift, and for 9 timestep pairs. For the 
SIS calculation the dynamical friction time gradually be- 
comes more of an overestimate as the simulation progresses 
in time. At high redshift z ^ 1, more than 65% of all subha- 
los spiral into their hosts in a time comparable to or longer 
than their calculated SIS dynamical friction time. At low 
redshifts, the SIS dynamical friction is typically an overes- 
timate, sometimes by more than two orders of magnitude. 
The use of the best-fit NFW profile brings about a somewhat 
worse agreement at high redshift. This is due to the lower 
halo concentration at high redshift that causes subhalos to 
be classified as "cD", residing within the R s of their hosts. 
The definition of r(tf) in such cases is more ambiguous and 
may not reflect the exact collision radius. Moreover, at high 
redshift a "cD" case may occur for progenitors of rather 
similar masses where the dynamical friction approximation 
breaks down anyway. In the SIS calculation for the "cD" 
cases, the steeper density profile near the center (oc r~ 2 ) 
tends to shorten the dynamical friction time estimate and 
partly compensates for the method inaccuracies. However 
at low redshift, where such degeneracies are rare, the NFW 
estimate does better than the SIS estimate and by z ~ 
it overestimates the time for only ~ 50% of all cases, in 
comparison to the ~ 70% in the SIS case. 

This spread of Tdf/At should partly be attributed to 
the fit errors, as at high redshift the halos tend to be less 
concentrated (BKS+) and the calculated dynamical friction 
time is affected at large radii by the modeled R s . The bigger 
errors for R s at high redshift propagate thus to an error in 
the dynamical friction time estimate. 

Semi-analytic models that use dynamical time estimates 
to predict mergers must take the errors in these estimates 
into account or else they underestimate the merger and col- 
lision rate both in absolute terms and more severely at low 
redshift relative to high redshift. 

Figure [Hi] summarizes Fig. |l8| and depicts the fraction 
of mergers between subhalos and hosts that occurred faster 
than the dynamical friction time pr ediction as a function of 
Klypin (1999] ) have analyzed the dy- 
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redshift. Kravtsov 



namical friction time-scale for clusters in these simulations 
using a different methodology, and found that at high red- 
shift the dynamical friction time-scale truly represents the 
cluster evolution whereas at lower redshift for a substantial 
fraction of the subhalos the dynamical friction timescale is 
longer than the actual merging timescale. 

Some of the mismatch between the Tdf and At may 
arise due to the invalidity of the circular orbit assumption. If 
most collision progenitors had initially radial orbits then the 
collision time should be smaller. The redshift dependence of 
the discrepancy can then be explained in terms of bigger 
infall velocities at lower redshift and domination of radial 
orbits near big groups and clusters. 

Figure ^ (top) shows the radial velocity of the subhalo 
progenitor normalized to the modeled circular SIS velocity 
(or equivalently v v [ r of NFW) at the lower redshift. The 
radial velocity should correlate with deviations from the dy- 
namical friction time predictions, especially the negative ra- 
dial velocities. Indeed, about 40% of the mergers for which 
the dynamical friction estimate failed to yield the right time 
scale, show radial velocities \v T \ > v v i T /s/2, whereas this frac- 
tion is only ~ 7% for mergers of correct time scale predic- 
tion. Interestingly, there is a marginal indication that even 
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Figure 19. Fraction of collisions where dynamical friction 
timescale overestimates the collision time as a function of look 
back time. 



positive radial velocities correlate with the breakup of the 
approximation. This may be interpreted as orbit eccentric- 
ity being an important reason for the deviation. At higher 
redshift (Fig. |20|, bottom) progenitors have more circular 
orbits (only 12% and 20% of \v t \ > iw/%/2 for correct and 
deviating mergers respectively) and therefore the deviation 
decreases. The difference between the time estimates by the 
SIS and NFW models is gradually amplified at low redshifts 
due to the larger concentration there (BKS+); the deviation 
of the simplified SIS profile from the realistic NFW profile 
grows with c v ir, which makes the SIS a weak approximation 
for high-c v i r halos. 

Figure |^ indicates that the deviations from the dynam- 
ical friction time estimate are correlated with the location of 
the progenitor subhalo (as identified in the output time just 
preceding the merger). Subhalos from the outskirts of the 
host tend to merge on a timescale much shorter than the 
dynamical friction time estimate. The discrepancies stem 
from such progenitors because they either spend longer ra- 
dial oscillation time there, or because they have just entered 
the host halo and are not yet virialized. A counter effect 
comes from deflections by other subhalos in the same host 
halo, which tend to make the dynamical fricti on timescale a n 
underestimate of the true spiraling-in time (Merritt 1983). 



Apparently this effect is overruled by the effects working in 
the opposite direction. 

Another possible cause for the mismatch between es- 
timated and observed Tdf may be the breakdown of the 
M su b/M hoBt <§C 1 assumption used in the derivation of r d f. 
We checked for a correlation between the deviation from dy- 
namical friction time estimate and the subhalo - host mass 
ratio. For the SIS estimate there does not seem to be any 
correlation between these two quantities. However all the 
deviating systems host - subhalo at high redshift are sys- 
tems for which M SI1 b/Mi, ost > 0.1 (the opposite is not true). 
Again, at high redshift when the concentration is low, and 
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Figure 20. Radial velocity in units of the host virial velocity 
as a function of the deviation from NFW dynamical friction time 
estimate. The dotted lines mark the 2 -1 / 2 ratio. For collisions 
between 0.03 > z > (top) and 2.9 > z > 2.3 (bottom). 



the subhalo mass is distributed over large radius, the as- 
sumption about a point-like subhalo with a fixed mass that 
does not affect the background potential breaks down, and 
hence the dynamical friction formula ceases to be valid es- 
pecially for similar masses of host halos and the spiraling 
subhalo. This effect is somewhat stronger when the actual 
density profile (i.e., NFW) is taken into account. Despite the 
fact that in the dynamical friction framework these effects 
tend to weaken the dynamical friction and thus make the dy- 
namical friction time estimate longer, other processes, e.g., 
two body relaxation are more common for "puffed" halos 
and may work to shorten the spiraling-in time scale. Note 
that by moving from the modeled M su h to M S ub(Rt), Tdf 
becomes longer (M sub <; M Bub (R t )) and thus the choice of 
using M su b cannot be the reason for the deviation. 
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Figure 21. Location of the subhalo progenitor (in units of the 
host i? v i r ) as function of the deviation from dynamical friction 
time estimate (SIS). 



7 SUMMARY 

We investigated the interaction rate of subhalos, which 
should be an important ingredient in the process of galaxy 
formation. This study is needed because virialized systems of 
objects evolve differently than unvirialized systems in which 
the objects are only weakly correlated. Also, the observed 
large fraction of galaxies in groups suggests that virialized 
systems of two galaxies or more are very common, and such 
high-density environments tend to be preferentially selected 
for observation. We thus computed the collision and merger 
rates within virialized systems in N-body simulations, and 
used the results to evaluate useful semi-analytic approxima- 
tions. The large sample of halos of a range of masses and 
multiplicities enabled a comprehensive statistical analysis, 
and the especially designed halo finder and collision identi- 
fier allowed a proper quantitative study of substructure evo- 
lution. Given the resolution limits of our simulations and 
halo finder, we could account for only the two most massive 
members of a group like the Local Group, but since these 
two members commonly dominate the group dynamics, the 
mass resolution is adequate for our purpose. 

We found that the comoving collision rate of subhalos 
is a substantial fraction to the total collision rate of halos 
in the ACDM simulation. At z ~ 4 the fraction of collisions 
that involve substructure is ~ 0.3, and this fraction grows 
to ~ 0.75 at low redshifts (Fig. |^ (top)). In physical coordi- 
nates, this corresponds to a collision rate oc (1 + z) 3 ~ 4 (Fig. 
Q). The evolution of the subhalo collision rate follows that 
of the total collision rate, namely a sharp rise between z = 4 
and 2, followed by a decline between z — 1 and 0.4. (Fig. ||). 
Even though the exact evolution is affected by the specific 
cosmology simulated and the results are valid for a specific 
mass range limited by the resolution of the simulation, the 
general behavior of a rise followed by a decline should be 
robust. We find that the collision rate of subhalos lags be- 
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hind the total rate by about 1 — 2 Gyr. It also exhibits a 
different behavior at low redshift, when many large virial- 
ized systems exist and contribute a significant fraction of 
the collision rate and a rise in the merger rate. 

The mass function of progenitors that participate in 
subhalo mergers is somewhat different from that of distinct 
halo collisions (Fig. |ll|), amplifying small differences in the 
tails of the mass functions for subhalos and distinct halos. 
The fraction of "major" interactions among those which in- 
volve at least one subhalo is only ^ 0.5, unlike the ~ 0.6 
fraction of major interactions among distinct halos. The 
moral is that accurate theoretical predictions of interaction 
rates should not be solely based on the halo mass function 
as derived from approximations such as Press-Schechter or 
Extended-PS, which ignore substructure. 

We found that the Makino-Hut formula can serve as 
a reasonably good approximation for subhalo merger rates, 
as long as host halos more massive than ~ 6 x 10 11 ft -1 Mq 
are considered. At lower masses, the MH formula overes- 
timates the merger rate by an order of magnitude. At high 
redshift, where less massive halos are abundant, the MH for- 
mula ceases to be a good description for the merger rate of 
the subhalo population. 

Predictions of the subhalo collision rate from kinetic 
theory are found to recover the actual rate for host halos 
more massive than ~ 10 13 ' 5 H'^Mq, but fail to do so, by 
one or two orders of magnitude, for less massive host halos, 
depending on their multiplicity. We presented a correction 
formula to match kinetic theory predictions with the simu- 
lated halos, using either the mass of the host, or the number 
of subhalos it contains. 

For interactions between subhalos and their hosts, we 
examined the accuracy of the dynamical friction timescale as 
estimated for circular orbits. In general, the estimate based 
on the NFW density profile is acceptable for about half 
the cases. It significantly overestimates the collision time in 
about 40% of the cases, and by more than an order of magni- 
tude in ~ 20% of the cases. One reason for this failure is that, 
at least at early times, the halos tend to be of low concen- 
tration and the high-mass host halos are of low abundance, 
making the assumption of M su b/Mh os t <§C 1 invalid in many 
cases. The deviation from circular orbits is clearly another 
reason for the failure of the dynamical-friction approxima- 
tion, as demonstrated at the most recent output time of the 
simulation. We find that, by some coincidence, the estimate 
based on the approximation of halos as isothermal spheres 
is acceptable at high redshift. This is good news for semi- 
analytic models that commonly use this approximation. At 
low redshifts, on the other hand, the concentration tends to 
be high, and the SIS predictions become poor. The moral is 
that, in order to mimic substructure dynamics in a cosmo- 
logical context, one needs to appeal to accurate simulations, 
or use approximations that were derived from such simula- 
tions. We have reasons to assume that these approximations, 
once expressed properly in terms of the host-halo mass, can 
be robust and not too sensitive to the particular cosmology. 
Proper collision rates can then be incorporated into galaxy 
formation scenarios and connect them to the evolution of 
large-scale structure. 
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